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1. INTRODUCTION 



Relatively little is known about queues for which 
interarrival times and/or service times are not independent; 
some studies are available however (cf. Jinlar [1967], 

Loynes [1962], Pearce [1967], Purdue [ 19 75 ]) . 

One reason for this seems to be the lack of tractable models 
for dependent sequences of random variables. Recently, models 
have been developed for sequences of dependent exponential 
random variables (cf. Jacobs and Lewis [1977]). These models 
are parametrically relatively simple and are in general not 

Markovian. Fortunately, it seems to be easy to model various 
types of dependence in queues using these sequences. 

The dependent sequence of exponential random variables 

is defined as follows. Let {e^} tie a sequence of independent 

random variables each with an exponential distribution with 

mean A ^ , 0 < A < 00 . Let { J n } and { } be independent 

sequences of independent {0,l}-random variables such that 

P{J = 1} = 1-8 and P{K =1} = 1-p where 0 < 8 < 1 and 
n n — — 

0 <_ p < 1 are fixed constants. For n = 1,2,... put 

( 1 . 1 ) 

where 

( 1 . 2 ) 

The sequence { } is called an EARMA(1,1) process (exponential 
mixed moving average autoregessive both of order 1) and ^ n _q 
is called the autoregressive part of X . 



X = 8 e 
n 



n 



+ J A . 
n n-1 



A = 
n 



pA n + 
n-1 



K e 

n n 



1 



We will always assume that has an exponential 

distribution with mean A ^ independent of all the other random 
variables. Under this assumption, Jacobs and Lewis [1977] have 
shown that {X^} is a stationary sequence with exponential 
marginal distribution having mean A ^ . Fuf thermore , the X r ' s 
are positively correlated with 

corr (X n ,X R+k ) = p k_1 (1-6) [ B ( 1 — p) + (1-6) p] , k > 1 . 

If 6 = 1, or 6 = 0 and p = 0, then {X n > is a 
sequence of independent random variables. The process {X^} 
is in general not Markovian although it is if 6 = 0 , in which 
case it is called an EAR(l) process (exponential autoregressive 
of order 1) . 

In this paper we will consider the simple case of a 
closed queueing network with two servers and a fixed number of 
cycling customers N. It is described as follows: when a customer 
finishes service at server 1 he joins the end of the queue at 
server 2; when he finishes service at server 2, he rejoins the 
end of the queue at server 1; the service discipline is first- 
in-first-out . 

This closed queueing network has been used in computer 
studies to model mul tiprogrammed computer systems (cf. Gaver 
and Shedler [1971]) and one is interested in obtaining, for 
example, the long run proportion of time one of the servers 
is idle, the average expected busy period of one server, and 
the average time it takes a customer to complete one cycle 



2 



of the network. There is some indication (cf. Lewis and Shedler 
[1973]) that the service times for one of the servers should 
be correlated. One object of this paper is to show that 
correlation of the service times does make a difference in the 
limiting behavior of the network. 

Let (respectively X n ) denote the nth service 

time for server 2 (respectively server 1) . There are many ways 
in which to use EARMA (1,1) processes to model dependence within 
each sequence of service times and cross correlation between 
the two sequences. In this paper we will assume that {S^} 
is a sequence of independent random variables each with an 
exponential distribution with mean y ^ , 0 < y < 00 and (X^ } 
is an EARMA (1,1) process independent of {S^} with mean X ^ , 

0 < X < 00 , and parameters 3 and p. Other queueing models 
using EARMA processes will be considered elsewhere. 

Let Z ( t ) denote the number of customers both waiting 
and being served by server 1 at time t. The process 
Z = {Z(t); t >_ 0 } takes the values {0,1,..., N). Let W(t) 
be the virtual waiting time at server 1 at time t; that is, 

W(t) is the sum of the service times of the customers in queue 
at server 1 at time t and the remaining service time of the 
customer currently being served. 

In the next two sections we will obtain limiting results 
for Z(t) and W(t) as t ->■ 00 for the case N = 2. We will 
show that the limiting distribution of Z (t) in the EARMA 
case is the same as in the case in which the service times of 
server 1 are independent; this result is also true for the case 
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N = 1 and seems to be due to the constraints imposed on the 
network by having only one or two customers. The long run 
average virtual waiting time is then computed for the case 
N = 2 and is found to be different from the independent case. 
This is because the positive correlation between the service 
times of server 1 tends to increase the virtual waiting time. 

In Section 4 we analyse the network for N > 2 customers 
and show the existence of a limiting distribution for Z (t) 
as t -*• °°. We then present some simulation results for the 
limiting distribution of Z(t) as t -*■ 00 for the case N = 5. 
The results show that the correlated service times do make a 
considerable difference in the limiting distribution. 



2 . THE TWO CUSTOMER CYCLIC QUEUE 

In this section we will study the process Z = {Z(t);t 0} 

and obtain the limiting distribution of Z (t) as t -*■ 00 for 
the queueing network of Section 1 with N = 2 . 

2 . 1 Preliminaries 

We define here three sequences which will be needed 
in the analysis. For concreteness we will always assume that 
Z ( 0 ) = 0. Let {T n } b e the increasing sequence of arrival and 
departure times for server 1 that includes all departure times 
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and those arrival times at server 1 that occur when server 1 

is idle. Let be the number of customers both waiting 

and being served just after time T . More precisely put 

Y n = Z ( 0 ) = 0 and T n = 0 and define Y , T , n > 1, 

0 0 n n — 

recursively as follows. If Y =0, let 

n 

(2.1) T ,, = inf { t > T : Z (t ) = 1}, Y = Z(T ,.+0) =1. 

n+1 n n+1 n+1 

If Y >0, let 
n 

(2.2) T n+1 = inf { t > T n :Z(t) < Z(t-)} , Y n+1 = Z(T n+1 +0) . 



For example, since Z(0) = 0 , is the time of the first 

arrival of a customer to server 1 and Y^ is the number of 

customers at server 1 just after time which must be 1; 

is the time of completion of the first service for server 1 

and Y 2 is the number of customers both waiting and being 

served at server 1 just after time ; the number could be 

0 or 1. In general, if Y > 0, then T is the time of 

completion of the next service for server 1 after time 

and Y n+ i -*- s num ^ )er customers both waiting and being 

served at server 1 just after time T . If Y =0, then 

J n+1 n 

T , 1 is the time of the next arrival to server 1 and Y , . = 1 
n+1 n+1 

Note that since we are considering the case in which N = 2, 

Y can only take the values 0 and 1. 
n 

Let denote the autoregressive part of the next 

service to be completed after time T . More precisely, put 

L rt = 0 and recursively define L, = inf{n > L, :Y =0}; 

U K+l K n 
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0. Then since 



L, is the index of the kth time that Y = 
k n 

Z (0 ) = °' A 0 = V A 1 = A q and A n = for 

L (k) < ri < L (k+1 ) . 

Since the service times at server 1 are dependent random 

variables, the dependence of { Z ( t ) ; t > T } on { Z ( t ) ; t < T } 

n — n 

is not only through Y n but also through the service times that 

were completed before time T . However, by (1.1) and (1.2), 

(Z(t); t > T^} is conditionally independent of {Z(t); t <_ T^} 

given (A ,Y ) . Further { (A ,Y ) } is a discrete time Markov 
^ n n n n 

process with state space (3R + x {0,1}, a (g x {0,1})) where 
1R + = [0,-] , R + denotes the Borel subsets of 3R + and 

a (R + x {0,1}) denotes the product a-algebra generated by R + 
and the subsets of { 0 , 1 }. 

Let 

P. . ( x , B) = P { A . € B , Y . = j j A = x, Y = i} 

13 n +1 n +1 J 1 n n 



and P(x,B) be the matrix whose i,j entry is P^(x,B) 

for i, j € {0,1}, x £ 3R and B € R + . The transition 

probabilities, P^j(x,B) are easily derived. If Y n = 0, 

then T . is the time of the next arrival after time T 
n +1 n 

and hence Y , . = 1 and A . . = A . If Y =1 and there is 
n+i n +1 n n 

an arrival during the service time “ T n » then ~ 1 / 

if there is no arrival during the service time “ T » 

then Y =0. Since {X } is an EARMA(1,1) process and 
{S n } is a sequence of independent exponential random variables 
with mean jj ^ we have the following matrix 
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0 




(2.3) P(x,B) 



-pX 

E[e X ;A 1 6B|A () =x] E[l-e 



-pX 



1 



;A 1 6B|A () = x] 



where <5 (B) is 1 if x € B and 0 otherwise. 

We will now define other processes related to (A ,Y ,T } 



which will also be used in the analysis. For T < t < T , , 
let Y(t) = Y , A ( t ) = A , and U(t) = t-T . Note that Z(t)=0 

if and only if Y(t) = 0; Z (t) =1 if and only if Y(t) = 1 

and there is no arrival at server 1 in the time interval 

( t-U ( t ) , t ] ; and Z(t) =2 if and only if Y(t) = 1 and there 

is an arrival in the interval (t-U(t),t]. Hence the limiting 
behavior of Z(t) as t -*■ 00 is related to that of 
(A(t) ,Y(t) ,U(t) ) . Furthermore, the limiting behavior of 
(A ( t) , Y (t) , U ( t) ) depends on that of (A n> Y n ). As a result, 
we will first compute the limiting distribution of (A n ,Y n ); 
then use it to compute the limiting distribution of 
(A ( t) , Y ( t ) , U ( t ) ) ; and finally compute the limiting distribution 
of Z (t) as t -*■ 00 . 

2.2 Limiting Properties of (A , Y ). 

Fix a Borel subset B of ]R with positive Lebesgue 
measure. From (1.2) it follows that {A } is a discrete time 



Markov process. To show that the process is recurrent note that 
expression (6.2) of Jacobs and Lewis [1977] for the kth order 



n n n 



n 



transition probability Q of {A^} implies that for each 
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dx there exists a k such that 



b > 0 and 0 £ 6 < / Ae 

B 

inf Q^(x,B) >6 . 

x € [0,b] 

Hence 

oo 

P( U {A € B} I A. = x) > 6 , x € 3R 

, m 1 0 — + 

m=l 

Therefore, by Proposition (5.1) of Orey [1971] 

P(A^ € B infinitely oftenlA^ = x) = 1 

for all x and thus {A^} is recurrent with respect to 
Lebesgue measure in the sense of Orey [1971, page 4], 

Since the service times for server 2 are independent 
with common exponential distribution, a similar argument shows 
that 

(2.4) P{A^ ^ B ' Y n = ^ infinitely often | Aq=x, Yg = i} = 1 

for i, j 6 {0,1}, x € 3R . Thus, by Theorem (7.1) of Orey 
[1971] there exists a possibly o-finite invariant measure 
{iTj(B); j = 0,1, B € R + } for the transition probability of 

{( vv>. 

The invariant measure satisfies the following system 
of equations 
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oo 



oo 



/ TT q ( dX ) P Q g ( X , B ) + / 7T 1 (dx) P 1 g(x,B) = 7T g ( B ) 

OO OO 

/ Tig (dx) Pg 1 (x,B) + / (dx) Pj^ ^XjB) = 77^ ( B) . 

Rewriting the two equations using the matrix (2.3) yields 

00 -yx 

(2.5) / 77^ (dx) E [e ; A 1 € B|A q =x] = 7T Q ( B) 



00 -yx 

(2.6) 7 Tq(B) + / 7i 1 (dx) E [ 1 - e 1 ;A 1 € B|A q =x] = (B) . 

Substituting the expression for 77 g (B) from (2.5) 
into (2.6), equation (2.6) becomes 



(2.7) / tt (dx) P(A C B | Ag=x) = tt (B) . 



By the result after (6.2) in Jacobs and Lewis [1977], equation 

(2.7) implies that 

(2.8) tt (B) = c / Xe _Ay dy 

1 B 

for any non-negative constant c. Substituting the expression 
for 77 ^ into (2.5) we have 
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(2.9) 



OO 



tt (B) = c / Ae" Ax E[e " 1 ; A € B|A =x] dx . 

0 l u 

We now want to choose c so that it is a probability and 

hence tt . (B) = lim ^ P{A € B, Y =i}. To this end we set 
l n -* 00 n n 

1 = tt (]R + ) + tt q (IR + ) 

00 -yx 

= c [1 + / Ae x E [e |A =x] dx] 

0 U 

= c [ 1 + A ( A+y ) _1 ] 

since, if A^ has an exponential distribution with mean A ^ , 
so does X^ . It now follows that c = ( A + y ) ( 2 A+i_i) ^ . 

This result will be used in the next subsection to compute the 
limiting distribution of ( A ( t ) , Y ( t ) , U ( t) ) . 

2 . 3 Limiting Properties of an Imbedded Semi-Markov Process 

Since ^ A n ,Y n ^ ,T n^ a s a Markov renewal process in the 
sense of (^inlar [1975], {[A(t),Y(t),U(t));t >_ 0} is a Markov 

process. Hence, { ( A ( t) , Y ( t ) ) ; t _> 0) is a semi-Markov process 
of the second type in the sense of Jacod [1973] . We will use 
the results of Jacod [1973] to compute the limiting distribution 
of (A ( t ) , Y ( t ) , U ( t) ) as t 00 which will then be used to 
compute the limiting distribution of Z (t) . 

Since {e } and {S^} are an< ^ e P en< ^ ent sequences of 
independent exponential variables, the process (A(t),Y(t)} 
is right continuous with left hand limits satisfying the 
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hypothesis R-3 on page 85 of Jacod [1973 ] concerning the set 
of discontinuity points of (A(t),Y(t)). Furthermore, 

{ (A ( t) , Y ( t) ) ; t >_ 0) is recurrent in the sense that, if 
B € R + has positive Lebesgue measure, then 



P{ / i, oy (s) 
0 { ^ } 



l B °A(s)ds = 



A_ = x , 



= i} = 1 



for i, j € {0,1} and x € IR where l 0 (x) is 1 if x € B 

T B 

and 0 otherwise. Therefore, Theorem III-10 on page 103 of 
Jacod [1973] applies to show that 



(2.10) v n (BxC) = d / ir n (dx) / e ys 1 (x) 1 (s)ds 
u 0 0 b c 



OO 00 

(2.11) v, (B x C) = d / tt (dx) / P{X > s|A =x} 1 (x) 1 (s)ds, 

i 0 1 o 1 U ^ 

B, C € R + , is an invariant measure for the transition function 
of the Markov process { (A (t ) , Y ( t) , U ( t) ) ; t >_ 0} for any 
d > 0. 

We will now choose d so that v^, i = 0,1, is a prob- 
ability. Put 

1 = VqUR^) + v 1 (3R^) = d [y -1 ir 0 ( 3R + ) + cA _1 ] 

= dc[A[y(A+y)] ^ + A ■*■] . 

Hence , 
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( 2 . 12 ) 



dc = Ay(A+y) (A 2 +Ay+y 2 ) ^ 



and d = Ay ( A+2y ) ( A 2 +Ay+y 2 ) 1 



2.4 The Limiting Distribution of Z (t) as t -*■ 00 

Using (2.10), (2.11), and (2.12) we can now compute the 

limiting distribution of Z (t) , the number of customers both 
waiting and being served at queue 1, as t -*■ °°. Recall 
that we are considering the case N = 2 and, therefore, Z(t) 
has the values 0, 1, or 2. 



From the argument after (2.3) it follows that 



(2.13) lim P{Z (t) =0 | Z (0) =0} = lim P{ Y (t ) =0 | Y ( 0 ) =0 } = v Q (HR 2 ) 




2 2 2 -1 
K ( A^ + Ay + v] 



(2.14) lim P{Z (t)=l | Z (0 ) =0 } 



lim E [e yU (t) ; Y ( t) =1 | Y (0 ) =0 ] 



t 00 



OO 



OO 



/ / v. (dx,ds) e ys 

0 0 



Ay [ A 2 + Ay + y 2 ] -1 



by (2.11) and (2.12); and finally 
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(2.15) lim P{Z (t) =2 | Z (0) =0} = lim E[(l-e yU ^ * ) ; Y ( t ) =l| Y ( 0 ) =0] 

^ oo ^ 00 

2 2 2 -1 
= y z [A z + Ay + y z ] . 

Note that the limiting distribution for the number of 
customers waiting or being served by server 1 with EARMA (1,1) 
service times is the same as if the service times for server 1 
were independent random variables each with an exponential 
distribution with mean A ^ (cf. Gaver and Thompson [1973]) . 

We feel that this is due to the constraints imposed on the net- 
work by having only two customers. This conjecture is shown 
to be true by the simulation results of Section 4 where it 
becomes clear that the result is not true if N > 2. 



3. THE VIRTUAL WAITING TIME FOR THE CYCLIC QUEUE WITH TWO 

CUSTOMERS 

In this section we will compute the long run virtual 
waiting time for server 1 in the case in which the service 
times for server 1 form an EARMA (1,1) process and there are two 
customers in the system. 

The virtual waiting time for server 1 at time t, W(t), 
is the sum of the service times of the customers in queue at 
server 1 at time t plus the remaining service time of the 
customer currently being served. To define W(t) more precisely, 
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let N(t) = sup{n:T n £ t) . If Z(t) is 1 or 2, then 
V(t) = T N ( t ) + i ~ t is the time from t until the completion 
of the current service. If Z(t) = 2, then T N(t)+2~ T N(t)+l 
is the length of the service time for the customer who arrived 
at server 1 in the time interval (t-U(t) ,t] . Hence, the 
virtual waiting time at time t is defined formally as 





1 0 


if 


Z(t) = 0 , 


(3.1) W(t) = * 


| v(t) 


if 


Z(t) = 1, 




V(t) + T N ( t ) +2 - 


T if 

N(t)+1 


Z(t) = 2 . 


From (3.1) 


it follows that 


{ (A ( t ) , Y ( t ) ,Z(t) ,U(t) ,V(t) 



t _> 0} is a Markov process. Further, the process is recurrent 
in the sense of Azema, Duflo, and Revuz [1969]. Thus by the 
ergodic theorem in Section (3.1) of that paper, there exists a 
constant W such that 



(3.2) W = lim i / W (s) ds 

t -► 00 0 

P(-|A(0)=x, Z (0 ) =0 ) -almost surely and 

(3.3) W = lim ± E[/ W ( s ) ds| A ( 0 ) =x , Z(0)=0] 

t - °° t 0 

for almost all x. We will use (3.3) to compute W. 

First, from (3.1) and the argument after (2.3) 
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(3.4) E[W(s) |A(0)=x, Z(0)=0] 

= E[V(s)e" yU(s) ; Y (s) =1 | A (0) =x, Z(0)=0] 

+E[ (V(s) + (T n(s)+2 - T N(s)+1 ))(l-e _MU(s) ) ;Y(s)=l|A(0)=x,Z(0)=0] 
= E[V(s) ;Y(s)=l |A(0)=x; Z(0)=0] 

+E[(T N(s)+ 2 " T N(s)+1 )(l-e _pU(s) );Y(s)=l|A(0)=x f Z(0)=0] 

We will first compute 
1 t 

(3.5) lim - / E [ V ( s ) ; Y(s)=l|A(0) = x, Y(0) = 0] ds . 

t -*■ 00 0 

Since (A ( t ) , Y ( t ) , U ( t ) ) has a limiting distribution as t °°, 
by the proof of Theorem III-10 on page 105 of Jacod [1973] , 
we have that 

p .(D) = lim P{ Y ( t ) =1 , (A(t) ,U(t) ,V(t) ) € D | A ( 0 ) =x , Y ( 0 ) =0 } 

1 t 00 

exists and equals 



CO oo 





cd 


/ 


Ae Az dz / P{X 


€ dy|A =z} / 


dx 1 (z,x,y-x) 






0 


0 


0 




where 


D 


is 


a Borel subset of 


3 

3R + . By the 


ergodic theorem of 


Azema , 


Duf lo 


and Revuz , (3.5) 


equals 




(3.6) 




/ 


(dz , dy , ds ) s = 


p (A + p) [A (A 2 + Ap + p 2 ) ] 1 . 



By (1.1) and (1.2), the second term of (3.4) 



15 



(3.7) E[(T N(s)+2- T N(s)+l ))(1_e ^* U(S) ^ ;Y(s)=l|A(0) =x , Z(0)=0] 



= e[3£ n(s)+2 (1_e yU(s) ); Y ( s ) =1 | A (0 ) =x , Z (0) =0 ] 



+ (l-3)E[pA(s) (1-e yU(s) ); Y(s)=l|A(0)=x, Z (0 ) =0 ] 



+ ( 1 — 3 ) (1-p) E [ £ N(s)+ i( 1 - e yU(3) ) ;Y(s)=l|A(0)=x,Z(0)=0] 



Since £ n( s j+2 i s independent of (Y(s),U(s)) 

(3.8) lim i/ E[£ (l-e" yU(s) ) ;Y(s)=l|A(0)=x,Z(0)=0]ds 

t -»• OO t 0 W 

-1 1 t 

= A lim i / P{Z(s)=2|A(0)=x,Z(0)=0}ds 

t -* 00 t 0 



2 2 2 -1 
- y [A (A Z + Ay+y ) ] 



by (2.15) and the previously cited ergodic theorem. 
The ergodic theorem also implies that 



(3.9) lim i- / E [A ( s ) ( 1-e yU ( S ) ) ; Y ( s ) =1 1 A ( 0 ) =x , Z (0 ) =0] ds 



t -* 00 0 



= / v (dz,dx)z(l-e yx ) 

2 1 

CO CO 

= cd / Ae ^ Z z dz / P{X > x|A =z) (1-e yX )dx 

0 0 1 u 
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oo oo 

= cd / Ae ^ Z z dz / [3P{Be, > x}+ (1-3) P{ Be, > x-z} ] (l-e -yx ) dx 

0 0 1 1 

= [X [ X 2 -t-Xu + n 2 ] (A+y) (A+By) ] _1 

xy 2 [3A 2 +2Ay+B (2Ay+2y 2 -A 2 ) + B 2 (-2 Ay-y 2 - X 2 ) ] 

by (2.11) , (1.1) and (1.2) . 

We now turn our attention to the computation of 

(3.10) lim i- f E[e N(s)+1 (l-e" pU(s) ) ; Y (s) =1 | A (0 ) =x, Z (0 ) =0] ds , 

First, note that (A ,Y ,e ,.) is a discrete time Markov process 

n n n+l 

with limiting distribution 77 ^(BxC) = tt^(B) / Xe ^ x dx for 
B, C 6 R + . Hence, by similar arguments to those in subsection 2.3, 
(A (t ) , Y ( t ) , , U ( t) ) has a limiting distribution v ^ . 

Further, for B, C, D € R , by Theorem III-10 of Jacod 

v. (BxCxd) = d / 77. (dz) /Xe ^ x dx / P (x > s|A =z, e =x)ds 
1 B C Dl 1 

= d fu, (dz) f Xe -Xx dx /{ B1 . ( 6x) + (1-B) 1 . , (Bx+z))ds 

B 1 C D S/ [ ' ’ 

for some constant d by (1.1) and (1.2). 

To evaluate d, note that 

(BxD) = V^B x 3R + x D) 

= dc / Ae~ >z dz / P{X > s|A =z) ds . 

B D 1 1 
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Hence, from (2.8) and (2.11), d = d. By the ergodic theorem 



again, (3.10) equals 



(3.11) / v (dz,dx,ds) x ( 1-e ps ) 

= (A(A 2 +Ap+p 2 ) (A+p6) ] _1 p 2 [A(l+6+6 2 ) + p(B+6 2 )] 



after some simplification. 

Putting together (3.3)-(3.11) we obtain after some 
simplification 



(3.12) [A (A 2 +Ap + p 2 ) ] _1 W 



= pA + 2p" 



+ (1-6) p 2 [ (A+p) (A+6p) ] 1 { [8A(A+p)+6 2 (A+p) 2 ] 

+p[2A(A+p)+6(-2A 2 +p 2 )-26 2 (A+p) 2 ] } 



Putting y = Ap ■*■ , the traffic intensity, (3.12) becomes 



(3.13) AW = (2+y) (y 2 +Y+1) 1 + (1-6) [ (Y 2 +Y+l) (Y+l) (Y+3) J" 1 

x { [ By (Y+l) +6 2 (y+1) 2 ]+p [2y (Y+l) +6 (-2 y 2 +D-26 2 (y + D 2 ] > 



If 6=1, then the service times for server 1 are 
independent exponential random variables with mean A ^ and 



(3.14) 



A W. nd = (2+y) (y 2 +Y+D 1 
18 



*-> aJ 

as expected in this case. Note that W >_ W i nc j and, hence, 
the positive correlation of the service times increases the 
average virtual waiting time. Also there is a such that 

/V 

for 3 < 8 q increasing p increases W; while for 8 > 6 Q 

/ v/ 

increasing p decreases W. The value Bq is a solution to 
the quadratic equation 

(3.15) -2(y+l) 2 B 2 + (-2y+l)B + 2y(y+l) = 0 

and hence depends on the traffic intensity y . 

To give an idea of the effect of different 8 and p 
on W we give Table 1 whose entries are differences between W 

r^J _ nj 

and W. , for various values of 8 and p for A=p = l(W. =1 
ind ind 

in this case) . 



p\B 


0 


.1 


. 3 


.5 


. 7 


.9 


0 


0 


.03 


.09 


.11 


.10 


.04 


.1 


.05 


.07 


.10 


.11 


.09 


.04 


.3 


.15 


.15 


.14 


.12 


.08 


.03 


.5 


.25 


.23 


.18 


.13 


00 

o 

• 


.03 


.7 


. 35 


. 30 


.21 


.13 


.07 


.02 


. 9 


.45 


.38 


.25 


.14 


.06 


> — 1 
o 

• 


TABLE 1. 


Values 


of W - 


"ind 


for the 


case 


A = y = 
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Note that it is the autoregressive part of the EARMA(1,1) 

/V# 

service time that causes the most change in W. As expected, 

/V 

for B small W increases with increasing p; for B large 
W decreases with increasing p; the value at which the change 
occurs, Bq , is .55 in this case. 



4. THE CYCLIC QUEUE WITH A FINITE NUMBER OF CUSTOMERS 

In this section we will study the queueing network of 

Section 1 in the case in which there are N >_ 3 customers. We 

will use without mention previous notation adapted to the 

present case. As before, we are interested in the limiting 

distribution as t -*■ 00 of Z(t), the number of customers either 

waiting or being served at server 1 at time t. Again, we first 

consider the discrete time Markov process {(A ,Y )} which now 

n n 

has state space (]R + x { 0 , 1 , . . . , N-l } , a ( g + x { 0 , 1 , . . . , N-l } ) ) 
where p(R + * { 0 , 1 , . . . , N-l } ) denotes the product a-algebra 
of the Borel subsets of 3R + and the subsets of { 0 , 1 , . . . , N-l } . 
Let P ± j (x , B) = P(A n+1 t B, Y n+1 =j|A n =x, Y n =i} for B € R + , 
x € ]R + , and i, j € { 0 , 1 , . . . , N-l } . If P(x,B) denotes the 
matrix whose (i,j) -entry is P„(x,B), then it is not hard to 
show that 
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0 


6 (B) 

X 


0 ... 


. . 0 


0 


b Q (x , B ) 


b 1 (x,B) 


b 2 (x , B ) 


2 ^ x ' B ^ 


a N _i(x,B) 


0 


b Q (x,B) 


h 1 (x,B) 


^N-3 ^ 


a N _ 2 ( x,B) 


0 


0 


b 0 (x,B) 


b N-4 ^ X ' B ^ 


3 ^ x ® ^ 



P (x , B ) = 



0 0 0 



b Q (x , B) a 1 (x,B) 



where 



b k (x , B) 




(MX ) k 

ITT A i e B IV X 



and 



a k (x ' B) = l b n ( x / B ) • 

n=k 



By similar methods to those of subsection 2.2 one can show 
that there exists a possibly o-finite invariant measure for P 
which satisfies the following system of equations for B £ R + 



(4.1) 



(4.2) 



/ TT, (dx)b (x,B) =TT (B) 
0 U 



00 



oo 



00 



/ 7T n (dx) 6 (B )+/ tt (dx)b (x,B)+/ tt (dx)b (x,B) =tt. (B) 

0 Xq 1 0 



21 



oo 



00 



(4-3) 



/ tt, (dx) b (x, B) + / 
0 1 0 



tt_ (dx) b (x, B) + / 

z 1 o 



ir 3 (dx) b Q (x ,B) =tt 2 (B) 



(4.4) / tt. (dx) a , (x,B)+/tt 7 (dx) a (x,B) + - ♦ • + / 

0 x 0 0 



7T N _i(dx)ai(x,B)=^ N _i 



(B) 



Substituting the expression for tTq(B) °f (4.1) into (4.2) 
and then adding all the equations except (4.1) together yields 
the equation 



/ tTT^+iT 2 +* • *+tt n _^] (dx) 



P{A € B | A q =x} 



[tt 1 +tt 2 + * • • +7T n _ 1 ] (B) . 



By the result after (6.2) of Jacobs and Lewis [1977] we have 



[tt^ +•••+ ir N _^] (B) = c ' / Xe dx 

B 

for any constant c'. Hence, we can and will choose c' so 
that the invariant measure it is a probability. 

We will now consider the limiting distribution of Z (t) 
as t -*■ 00 . Let K ( t ) be the number of service completions by 
server 2 in the time interval (t - U(t) ,t] . Since the service 
times of server 2 are independent each with the same exponential 
distribution { (A (t) , Y (t ) , U ( t) ,K (t ) ) ; t >_ 0 } is a Markov process. 
Note that (Z(t)=0) = {Y(t)=0}; 
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{ Z ( t ) =i } 



0 < i < N; 



and 



= U {Y(t)=k,K(t)=i-k} , 
k=l 

N-l 

{ Z ( t ) =N } = U {Y (t) =k,K(t) =N-k} . 

k=l 



Thus, using the techniques of subsections 2.3 and 2.4 one can 
show that 



v (0 ) = lim P{Z (t ) =0 | Z (0)=0,A(0)=0} = dy 1 tt (JR ) ; 

t -> °° 



v (N) = lim P { Z ( t ) =N | Z ( 0 ) =0 , A ( 0 ) =0 } 
t -*■ 00 



N-l 



= a l [ tt (dy) / P{X >s|A =y} l e 
k=l 0 K 0 1 1 j=N-k 



-ys (ys) ] 



-i ! 



ds ; 



and 



v (i ) = lim P{Z(t)=i|z(0)=0,A(0) =x } 



1 00 



d l / 77 k (dy) / P{X >s|A =y} e ys , 

k=l 0 0 Kl ' 



i-k 



ds 



for some positive constant d. 

It seems difficult to solve (4.1)-(4.4) for tt^, 
i = 0,...,N-1. Hence, we are unable to obtain explicit expressions 
for the limiting distribution of Z (t) as t -*• 00 as we could in 
the case N=2 . We will, however, give some simulation results 
to indicate the limiting behavior of Z (t) as t ->• 00 for higher 
values of N. 

The simulation is based on the following observation. 

Since the service times of server 2 are independent and exponentially 
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distributed { (A ( t ) , Y ( t ) , U ( t ) , Z ( t ) ) ; t ^ 0} is a Markov process 
that satisfies the hypotheses of the ergodic theorem of Az£ma, 
Duflo and Revuz [1969]. Hence, 

t 

lim — / lr . i °Z(s) ds = v(i) , i = 0,...,N 

t -> co t o x11 

almost surely P(*|A(0)=x, Z(0)=0). Further, if x M is the 

time of completion of the Mth service time for server 1, then, 

since t -*■ 00 as M -*■ m almost surely with respect to 
M 

P(. | A(0)=x,Z(0)=0) 

1 TM 

(5.5) lim / 1 r . |°Z (s)ds = v(i), i = 0 , . . . ,N 

M -* 00 M 0 ll/ 

almost surely with respect to the same probability. 

The following results are from a simulation designed 
by Professor P. A. W. Lewis at the Naval Postgraduate School, 
Monterey, California for a closed queueing network having 
EARMA (1,1) service times. The simulation we used consists of 
the computation of 

1 r X M 

(5.6) / lr.,oZ (s)ds 

t m 0 x 1 

for M equal to the 5000th service time. The computation was 
repeated for 100 independent replications and the sample mean 
and variance over the 100 replications were computed. All 
runs were performed on an IBM system 360/67 computer at the 
Naval Postgraduate School using the LLRANDOM package (Learmonth 
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and Lewis [1973]) which generates numbers according to the 
scheme given by Lewis, Goodman, and Miller [1969] and 
exponentially distributed random numbers using the Marsaglia 
" rectangle-wedge- tail " method. Tests of the random number 
generator are given in Learmonth and Lewis [1974]. 

Table 2 gives the results of a simulation of the queueing 
network for the case N=5 and A=y=l. The entries in the table 
are the sample means of (5.6) over the independent replications 
for different values of $ and p. The numbers in the 
parentheses are estimates of the standard deviations of the 
estimate of v(i). For comparison, note that the limiting 
distribution in the case in which server 1 has independent 
exponential service times and A=y=l is 

(5.7) lim P { Z ( t ) =i } = 1/6 = 0.1666, i = 0,...,5 , 

t 00 

(cf. Gaver and Thompson [1973]) . This corresponds to the case 
6 = 0, p = 0 in the first line of the table. 

Note that again it is the autogressive part of the 
service times that causes most of the change in the estimates 
for the limiting distribution of Z (t) as t ->• 00 . The 
positively correlated service times increase the probabilities 
of server 1 or server 2 being idle; they also increase the 
probability of all customers being in one or the other service 
center. This seems to be due to the fact that, if p is large 
and 3 is small, then having a large service time at one time 
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i 


0 


1 


2 


3 


4 


5 


3= 0 


p= 0 


. 17 ( .009) 


. 17( .007) 


. 17( . 005) 


.17 ( .006) 


.17 ( . 007) 


. 17 ( .009) 




.1 


. 18 ( . 01) 


. 16 ( .007) 


. 16 ( . 005) 


. 15 ( .005) 


. 17 ( .007) 


.17 ( .01) 




.5 


.21(.01) 


. 13( . 006) 


. 12 ( .006) 


. 12 ( . 005) 


. 18( .007) 


.21(.0l) 




.9 


. 26 ( .02) 


.ll(.Ol) 


.079 ( . 006) 


.088 (.006) 


.20( .01) 


. 26( .03) 


B-.l 


p= 0 


. 1 7 ( - 01) 


. 16( .008) 


. 16 ( .005) 


. 16 ( .005) 


.17 ( .007) 


.17 (.01) 




.1 


. 18 ( .01) 


.16 (.008) 


.16( .005) 


. 15( .005) 


. 17 ( .007) 


. 17 ( .01) 




.5 


.21(.01) 


. 15( .007) 


. 13 ( . 006) 


.13( .005) 


. 18 ( . 007) 


. 20( . 01) 




.9 


. 2 5 ( .02) 


. 13 ( .01) 


.092 ( .006) 


. 10( . 005) 


.18( .01) 


.25 (.03) 


3= . 5 


p= 0 


. 18 ( .01) 


. 17 ( .007) 


. 16 ( .005) 


.15 ( . 006) 


.16 ( .007) 


. 18( .01) 




.1 


. 18 ( . 01) 


.17 ( .007) 


.16( .005) 


. 15 ( .006) 


. 16 ( .007) 


. 18 ( .01) 




.5 


. 19 ( . 01) 


. 17 ( . 007) 


. 15 ( .006) 


. 15 ( . 006) 


.16( .007) 


.18 (.01) 




.9 


. 20 ( . 02) 


. 16 ( .01) 


. 14 ( .007) 


.14 (.006) 


.16( .01) 


. 20 ( .02) 


$=.7 


P= o 


. 18 ( .01) 


. 17 ( .007) 


. 16 ( . 005) 


. 15 ( .005) 


. 16 ( .007) 


.18( .01) 




.1 


. 18( . 01) 


. 17 ( .007) 


. 16 ( .005) 


. 15 ( . 006) 


. 16 ( . 007) 


. 18( .01) 




.5 


. 18 ( .01) 


.17 ( .008) 


. 16 ( .005) 


. 15 ( .006) 


.16 (.008) 


. 18( .01) 




.9 


. 18 ( .01) 


. 17 ( .009) 


. 16 ( . 006) 


. 16 ( .006) 


. 16 ( .009) 


. 17 ( .01) 


3=.9 


P= o 


. 17 ( .01) 


. 17 ( .007) 


.16 (.005) 


. 16 ( .006) 


. 16 ( . 007) 


. 17 ( .01) 




.1 


. 17( .01) 


.17 ( .007) 


.16 (.005) 


.16( .007) 


.16( .007) 


. 17( . 01) 




.5 


. 17 ( . 01) 


.17 ( .007) 


. 16 ( .005) 


.16( .006) 


.16(.008) 


. 17( .01) 




.9 


. 17 ( .01) 


. 17 ( .008) 


. 17( .005) 


,16( .006) 


.16 (.008) 


. 17 ( .01) 



TABLE 2. Estimates for the limiting distribution for the number of customers 

in queue for N=5 customers and X = 1 when server 1 has EARMA(1,1) 
service times and server 2 has independent exponential service times. 
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implies that the service times will tend to be large for a 
while; having a small service time implies that the service 
times will tend to be small for a while. The dependent service 
times appear to have little effect on the long run average 
number of customers that are waiting or being served at server 1; 
the average is 2.5 for all values of (B,p) except the values 
(0,.9), (.1,.5), and (.1,.9) for which it is 2.6. This 

invariance reflects the fact that the change in the limiting 
distribution (v(i)} due to the dependence is somewhat 
symmetric in i. 

Finally, we present in Table 3 the results of a simulation 
to investigate what happens if the service times of both servers 
are independent EARMA (1,1) processes. In that simulation 
the number of customers N is 5 and the parameters for the 
EARMA (1,1) service times for server 1 are B-^ = .1 and = .9. 

The parameters and P 2 for the EARMA (1,1) service times 

for server 2 are allowed to take on several values. The 
expression (5.6) is computed for M = 10,000 for 750 independent 
realizations. The entries in the table are the values for this 
simulation of the same quantities as in Table 2. 



i 


0 


1 


2 


3 


4 


5 


B 2= 


.1, 


O 

II 

CM 

Q. 


.249 (.02) 


.130( .006) 


.090( .004) 


. 100 ( . 004) 


. 182 ( .008) 


,249( .016) 


V 


.1, 


P 2 =. 5 


.265 ( .017) 


. 142 ( . 007) 


.079 ( .003) 


.086 ( .004) 


. 164 ( . 007) 


. 265 ( . 017) 


6 2= 


.1, 


P 2 =-9 


. 290 ( .021 ) 


.149 (.009) 


. 060 ( . 004) 


. 060 ( .004) 


. 150 ( . 009 ) 


.290( .021) 



TABLE 3. Estimates of the limiting distribution for the number 
of customers in queue for N=5 customers and A=p=l . 
The two sequences of service times are independent 
EARMA (1,1) processes. 
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Note that the positive correlation of both sequences of 
service times tends to make the probability of server 1 being 
idle and the probability of all customers being at server 1 
larger than when only one sequence of service times is 
positively correlated. The change in the limiting distribution 
v(i) is again somewhat symmetric in i. The long run average 
number of customers waiting or being served at server 1 is 
2.5 if 62 = *1/ P 2 = *9 and 2.6 for the other values of 
( 3 2 ' P 2 ^ * 

It is clear from both simulations that the limiting 
distribution for number of customers at server 1 is quite 
sensitive to serial correlation in the service times. The 
simulations indicate that perhaps v(0) = v(5) in the case 
A = y = 1 . It can in fact be shown that in general 
A [1 - v ( 0 ) ] = y [1 - v (N) ] . 

5. CONCLUSIONS 

In this paper we considered one scheme for using EARMA 
processes to model dependence in queues. We find that the 
introduction of dependence does affect the limiting behavior 
of the queue. There are, of course, many other schemes and 
some of these will be considered elsewhere. Two advantages 
of using the EARMA processes in queues are the ease of intro- 
ducing dependence in the queue and the ease of simulating 
the processes. The major drawback to using EARMA processes 
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in queues is the difficulty of obtaining the exact analytic 
results. It is expected, therefore, that approximation 
techniques and simulation will be of major importance in 
analysing these queues. 
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